Onset of a Quantum Phase Transition with a Trapped Ion Quantum Simulator 
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A quantum simulator is a well controlled quantum system that can simulate the behavior of another quan- 
tum system which may require exponentially large classical computing resources to understand otherwise. In 
the 1980s, Feynman proposed the use of quantum logic gates on a standard controllable quantum system to 
efficiently simulate the behavior of a model Hamiltonian (IIU- Recent experiments using trapped ions jsj-Q] 
and neutral atoms |7] have realized quantum simulation of Ising model in presence of external magnetic fields, 
and showed almost arbitrary control in generating non-trivial Ising coupling patterns |4]. Here we use laser- 
cooled trapped ^^^Yb"*' ions to simulate the emergence of magnetism in a system of interacting spins (sl- fioll 
by implementing a fully-connected non-uniform ferromagnetic Ising model in a transverse magnetic field. To 
link this quantum simulation to condensed matter physics, we measure scalable correlation functions and order 
parameters appropriate for the description of larger systems, such as various moments of the magnetization. 
By increasing the Ising coupling strengths compared with the external field, the crossover from paramagnetism 
to ferromagnetic order sharpens as the system is scaled up from A'^ = 2 to 9 trapped ion spins. This points 
toward the onset of a quantum phase transition that should become infinitely sharp as the system approaches 
the macroscopic scale. We compare the measured ground state order to theory, which may become intractable 
for non-unifoiTn Ising couplings as the number of spins grows beyond 20- 30 [ 1 1,1 and even NP complete for 
a fully-connected frustrated Ising model fl2ll . making this experiment an important benchmark for large-scale 
quantum simulation. 



We find the ground state of the following Hamiltonian for 
interacting spin— 1/2 systems : ifTsIl 

H^-^Y.As'^^i-BY.^l (1) 

i<j i 

where ct^ is the Pauli matrix for the i^^ spin (a = x,y, z and 
i = 1,2,..., N) and J^.j > is the ferromagnetic (FM) Ising 
coupling matrix, with J — (Jij). Our experiment is per- 
formed according to adiabatic quantum simulation protocol 
(14'] where the dimensionless coupling B/\J\ is tuned slowly 
enough so that the system follows instantaneous eigenstates 
of the changing Hamiltonian illl. As B/\J\ 00, the 
ground state has all spins polarized along the magnetic field, 
or paramagnetic along the Ising direction X. In the other limit 
B/\J\ — 0, the spins order according to the Ising couplings 
and the ground state is a superposition of FM states | tt ... t) 
and I II ... I) where | t) and | 4,) are eigenstates of ax- We 
characterize the magnetic order in the system by measuring 
various correlation functions between all N spins, including 
the probability of FM occupation and the second and fourth 
moments of the total magnetization. 

We represent each spin-i system by the hyperfine clock 
states ^S'l/alF = 0,mp = 0) and \F ^ l,m,p ^ 0) of 

i^iYb+ separated by i^hf = 12.642821 GHz, which are de- 
noted by the eigenstates and of cr^, respectively. These 
states are detected by standard spin-dependent resonant fluo- 
rescence on the cycling ^6*1/ 9 to '^Pi/2 transition at 369.5 nm 
using a photomultiplier tube 01511 . The ions are trapped along 
the z— axis of a three layer linear Paul trap (see Fig. [TJ) with 
center of mass (CM) vibrational frequencies of ly^ — 4.748, 



= 4.641, and = 1.002 MHz along the x.y (trans- 
verse) and z (axial) directions, respectively |[l6ll . The modes 
of motion along x are cooled to near their vibrational ground 
states and within the Lamb-Dicke regime. Off -resonant laser 
beams address the ions globally, driving stimulated Raman 
transitions between the spin states and also imparting spin- 
dependent forces exclusively in the x-direction, as depicted 
in Fig.[TJ) 10-0]. The Raman beams contain a "carrier" beat- 
note at frequency i^hf, which provides an effective uniform 
transverse magnetic field B. Raman beatnotes at frequencies 
vhf ^ are near motional sidebands and generate a spin- 
dependent force. The rf phase difference between the carrier 
beatnote and the average beatnote of the sidebands is set to 
7r/2 so that the magnetic field is transverse to the Ising cou- 
plings (Eq 1) m laf. We suppress direct sideband (phonon) 
excitation by ensuring that the beatnote detuning /i is suffi- 
ciently far from any normal mode frequency [4]. This requires 
that l/i — Vm\ » rji^m^i, where rji^ra is the Lamb-Dicke pa- 
rameter of the zth ion and mth normal mode at frequency v„i 
(with vi — Vx denoting the CM mode), and il^ = g^/ A is 
the carrier Rabi frequency on the zth ion. Here gi is the single 
photon Rabi frequency of i*'* ion and A is the detuning of the 
Raman beams from the — ^ P1/2 transition. 

This results in an Ising interaction between the spins with 
control parameter fi that dictates the form of the coupling ma- 
trix 1,4J 

N 

In the experiment A ^ 2.7 THz, Hi ^ 370 KHz and we ex- 
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Figure 1. Experimental schematics : a. Schematic of the three-layer linear rf (Paul) trap, with top and bottom layers carrying static potentials 
and the middle one carrying rf. b. Two Raman beams globally address the ^^^Yb^ ion chain, with their wavevector difference (Afc) along 
the transverse (x) direction of motion, generating the Ising couplings through a spin-dependent force. The same beams generate an effective 
transverse magnetic field by driving resonant hyperfine transitions. A CCD image showing a string of nine ions (not in present experimental 
condition) is superimposed. A photomultiplier tube (PMT) is used to detect spin-dependent fluorescence from the ion crystal, c. Outline of 
quantum simulation protocol. The spins are initially prepared in the ground state of —B'^- Oy, then the Hamiltonian|T]is tumed on with 
starting field Bo 3> \J\ followed by an adiabatic exponential ramping to the final value B, keeping the Ising couplings fixed. Finally the X — 
component of the spins are detected. 



pect Ji.j/N ^ 1 KHz for the beatnote detuning ji such that 
jjL — vi Arji^ifti, as shown in Fig.|2^ . This beatnote cor- 
responds to 63 KHz blue of the CM mode frequency for 2 
ions and 30 KHz for 9 ions, as the Lamb-Dicke parameter 
r]i „i ~ This maintains roughly the same level of vir- 

tual phonon excitation for any number of ions. The expected 
Ising coupling pattern for a uniformly illuminated ion chain is 
shown in Fig.|2j3 for = 9 ions and the couplings are dom- 
inated by uniform contribution of the CM mode. The non- 
uniformity in the Ising couplings arises from other vibrational 
modes and a < 1% variation in fi^ across the ion chain for 
^ 100 fim x5/im wide gaussian Raman beams used in the 
experiment. 

The experiment proceeds as follows (Fig.[l];). We initialize 
the spins to be aligned to the Y— direction of the Bloch sphere 
by optically pumping to | l^iz ■■■ iz) and then coherently ro- 
tating the spins through tt/2 about the Bloch X— axis with a 
caiTier Raman transition. Next we switch on the Hamiltonian 
H with effective magnetic field Bq ^ 5\J\ so that the spins 
are prepared predominantly in the ground state. Then we ex- 
ponentially ramp down the effective magnetic field with a time 
constant of 80 /is to a final value B, keeping the Ising cou- 
plings fixed. We finally measure the spins along the Ising (X) 



direction by coherently rotating the spins through tt/2 about 
the Bloch F— axis before fluorescence detection. We repeat 
the experiment ^ lOOOA^ times and generate a histogram of 
fluorescence counts and fit to a weighted sum of basis func- 
tions to obtain the probability distribution P{s) of the number 
of spins in state (| t)), where s — 0, 1, ...N, as described in 
the methods section. 

We can generate several magnetic order parameters of in- 
terest from the distribution P{s), showing transitions between 
different spin orders. One order parameter is the average ab- 
solute il8il magnetization per site along the Ising direction, 

N 

"^a: — Jj J2 " 2s\P{s). For a large system, this param- 

s=a 

eter shows a second order phase transition, or a discontinuity 
in its derivative with respect to B/\J\. On the other hand, the 
fourth-order moment of the magnetization or Binder cumu- 

N / N \ ^ 

lant (BC) 5 = E " 2s)*P(s)/ E " 2s)^P{s) 
j—ir-L s=o Vs=o / 

B191 12011 becomes a step function at the QPT and should there- 
fore be more sensitive to the phase transition. We illustrate 
this point by plotting the exact ground state order in the simple 
case of uniform Ising couplings for a moderately large sys- 
tem (N — 100) in Fig. [3^. In figure [Sji-d we present data 
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Figure 2. Motional modes and Ising couplings: Transverse flT'l vibrational modes are used in the experiment to generate Ising couplings 
according to Eq. |2] a. Raman sideband spectrum of vibrational normal modes al ong transverse a;— direction for nine ions, labeled by their 
index m. The two highest frequency modes at i/i (CM mode) and 1^2 ~ — i'l ("tilt" mode) occur at the same position independent of the 
number of ions. The dotted and the dashed lines show beatnote detunings of fj, x ui + 30 KHz and « i^i + 63 KHz used in the experiment 
for A'^ = 9 and N = 2 ions respectively. CaiTier transition, weak excitation of transverse— j/ and axial— 2 normal modes and higher order 
modes are faded (light grey) for clarity, b. Theoretical Ising coupling pattern (Eq. ^ for N — 9 ions and uniform Raman beams. The main 
contribution follows from the uniform CM mode, with inhomogeneities given by excitation through the other nearby modes (particularly the 
tilt mode). Here, Ji,i+r ~ l/r"'^^ (r > 1), as found out empirically. For larger detunings, the range of the interaction falls off even faster 
with distance, approaching the limit Jij ^ — j\^ for /i ^ lyi 



for these two order parameters as B/\J\ is varied in the adia- 
batic quantum simulation, with the order parameters properly 
scaled to account for trivial finite size effects, as described 
in the methods section. Fig. [3}3 shows scaled magnetization, 
rhx for = 2 to iV = 9 spins that depict the sharpening 
of the crossover curves from paramagnetic to ferromagnetic 
spin order with increasing system size. The linear time scale 
indicates the exponential ramping profile of the (logarithmic) 
B/\J\ scale. Fig.[3};-d compares the two extreme system sizes 
in the experiment, N = 2 and iV = 9 and clearly shows the 
increased steepness for larger system size. The scaled mag- 
netization fhx is suppressed by ^ 25% (Fig. [3j3,c) and the 
scaled BC g is suppressed by ^ 10% (Fig. [3]i) from unity 
at B/\J\ = 0, predominantly due to decoherence from off- 
resonant spontaneous emission and additional dephasing due 
to intensity fluctuations in Raman beams during the simula- 
tion. 

We compare the data shown in Fig.[3};-d to the theoretical 
evolution taking into account experimental imperfections and 
errors discussed below, including spontaneous emission to the 
spin states and states outside the Hilbert space, and additional 
decoherence. The evolution is calculated by averaging 10, 000 
quantum trajectories. This takes only one minute on a sin- 
gle computing node for N = 2 spins and approximately 7 
hours, on a single node, for N ^ 9 spins. Extrapolating from 
this calculation suggests that averaging 10, 000 trajectories for 
= 15 spins would require 24 hours on a 40 node cluster, 
indicating the inefficiency of classical computers to simulate 
even a small quantum system. 

A faithful quantum simulation requires an excellent under- 
standing of errors and their scaling, especially when the un- 
derlying problem is otherwise intractable. We characterize er- 



rors in the cuiTent simulation by plotting the observed param- 
eter P(FM)=P(0)H-P(N) for iV = 2 to iV = 9 spins in Fig. |4] 
and comparing with theory. The parameter P{FM) involves 
only two of the 2^ basis states and is therefore more sensi- 
tive to errors compared with the order parameters fhx and g. 
For instance, ai B /\J\ =0 in Fig.[3}5-d and Fig.|4^-d we find 
that fhx and g do not change appreciably with system size, but 
P(FM) degrades to ~ 0.55 for = 9 spins from ^ 0.9 for 
iV = 2. 

There are several primary sources of experimental error 
Diabaticity due to finite ramping speed and eiTor in initial- 
ization is estimated to suppress P(FM) by ~ 3% for N = 2 
to ^ 8% for A^ = 9. This also gives rise to oscillations seen 
in the data (Fig. ^-d and Fig. |4^-d). A major source of er- 
ror is the spontaneous emission from Raman beams which 
amounts to a ^ 10% spontaneous emission probability per 
spin in 1 ms. Spontaneous emission dephases and randomizes 
the spin state and loosely behaves like a 'spin temperature' 
in this system, though the spins do not fully equilibrate with 
the 'bath' and the total probability of spontaneous emission 
increases linearly during the quantum simulation. In addition, 
each spontaneous emission event populates other states out- 
side of the Hilbert space of each spin with a probability of 
1/3. Spontaneous emission errors grow with increasing sys- 
tem size, which also suppresses P(FM) order with increasing 
A^, as seen in Fig.|4^-d. We theoretically estimate the suppres- 
sion of P(FM) due to diabaticity and spontaneous emission to- 
gether by averaging over quantum trajectories to be ^ 7% for 
N — 2 spins and ^ 26% for A^ = 9 spins. Intensity fluctua- 
tions on the Raman beams during the simulation modulate the 
ac Stark shift on the spins and dephase the spin states, which 
causes additional diabaticity and degrades the final ferromag- 
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Figure 3. Experimental results of adiabatic quantum simulation : a. Theoretical values of order parameters rhx and g are plotted vs B/l J| 
for a moderately large system (TV = 100 spins) witli uniform Ising couplings, in the case of a perfectly adiabatic time evolution. The effective 
ground state manifold reduces to A'^ + 1 dimensions in the total spin basis when the Ising couplings are uniform. The scaled Binder cumulant 
g approaches a step function near the transition point B/\J\ = 1 unlike the scaled magnetization rhx, making it experimentally suitable to 
probe the transition point for relatively small systems, b. Scaled average absolute magnetization per site, rhx vs B/\J\ (and simulation time) is 
plotted for A'^ = 2 to TV = 9 spins. As B/\J\ is lowered, the spin ordering undergoes a crossover from a paramagnetic to ferromagnetic phase. 
The crossover curves sharpen as the system size is increased from TV = 2 to TV = 9, anticipating a QPT in the limit of infinite system size. 
The oscillations in the data arise due to imperfect initial state preparation and non-adiabaticity due to finite ramping time. The (unphysical) 
3D background is shown to guide eyes. c. Magnetization data for TV = 2 spins (circles) is contrasted with TV = 9 spins (diamonds), near 
B/\J\ — 1. The data deviate from unity at B/\J\ = by ~ 20%, predominantly due to decoherence from spontaneous emission in Raman 
transitions and additional dephasing from Raman beam intensity fluctuation, as discussed in the text. The theoretical time evolution curves 
(solid line for TV = 2 and dashed line for TV = 9 spins) are calculated by averaging over 10,000 quantum trajectories, as described in the 
Methods section, d. Scaled Binder cumulant (g) data (circles and diamonds) and time evolution theory curves (solid and dashed lines) are 
plotted for TV = 2 and TV = 9 spins. At B/\J\ =0 the data deviate by ~ 10% from unity, due to decoherence as mentioned before. 



netic order. When we introduce a theoretical dephasing rate of 
0.3 per ms per ion in the quantum trajectory computation the 
predicted suppression of P(FM) increases to ^ 9% for TV = 2 
and ~ 50% for N = 9. 

Imperfect spin detection efficiency contributes ~ 5 — 10% 
error in P(FM). Fluorescence histograms for P(0) and P(l) 
have a ^ 1% overlap (in detection time of 0.8 ms) due to off- 
resonant coupling of the spin states to the level. This 
prevents us from increasing detection beam power or photon 
collection time to separate the histograms. Detection error 
in the data include uncertainty in fitting the observed fluo- 
rescence histograms to determine P{s), intensity fluctuations 
and finite width of the detection beam. 

The role of phonons in the results of quantum simulation 



is investigated both experimentally and numerically. We per- 
form another set of experiments with ^ — vi =63 KHz for 
all TV, which amounts to reducing the phonon excitation as the 
number of spins TV is increased, since the Lamb-Dicke param- 
eter rji^rn ^ We do not note any appreciable difference 
(beyond the margin of experimental errors) in spin population 
with the results reported here. In the presence of the applied 
effective magnetic field, phonon modes are coherently popu- 
lated and generally exhibit spin-motional entanglement. How- 
ever, these phonons do not alter the spin ordering and hence 
preserve spin-spin correlation, even if the entanglement be- 
tween spin states is partly destroyed when measurement traces 
over phonon states. Thus a standard GHZ type witness opera- 
tor will not show entanglement beyonds a few spins. 
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Figure 4. Suppression of P(FM) with increasing number of spins : a-d. Ferromagnetic order P(FM)=P(0)+P(N) is plotted vs B/\J\ for 
= 2 to A'^ = 9 spins. The circles are experimental data and the lines are theoretical results including decoherence and imperfect initialization. 
As this quantity includes only two of 2^ basis states random spin-flips and other errors degrade it much faster than the magnetization and 
Binder cumulant. The representative detection error bars are shown on a few points for each A'^. The P(FM) reduces from ~ 0.9 to ~ 0.55 as 
the system size is increased from two to nine. The principle contribution to this degradation is decoherence, predominantly due to spontaneous 
emission from intermediate states in the Raman transition and additional dephasing primarily due to intensity fluctuations in Raman 

beams. Shown in d is an estimated breakdown of the suppression of P(FM) from various effects for A'^ = 9 spins. Non-adiabaticity due to 
finite ramping speed, spontaneous emission and additional dephasing due to fluctuating Raman beams suppress P(FM) by ~ 8%, ~ 18% and 
~ 24% respectively from unity (B/\J\ — > 0). 



In this experiment we have qualitatively observed sharpen- 
ing of crossover curves which indicates the onset of a quan- 
tum phase transition as the system size increases. This scheme 
can be scaled up to larger number of spins where it is possi- 
ble to quantitatively estimate finite size effects e.g., scaling 
in critical exponents near the phase transition point [21|. The 
primary challenges in experimenting with larger system sizes 
include the spontaneous emission as described above, the re- 
quirement of larger optical power to maintain same level of 
Ising couplings, and non-adiabatic effects due to shrinking 
gap between ground and first excited state ^23, ^] . One solu- 
tion is to implement a high power laser with a detuning far 
from the energy levels, which would minimize sponta- 
neous emission while maintaining the same level of Ising cou- 
plings. This would also allow versatility in varying the Ising 
interaction (together with the effective external field) during 
the simulation, as the differential A.C. stark shift between spin 
states is negligible for a sufficiently large detuning. The coher- 
ence time increases in the absence of spontaneous emission, 
allowing for a longer simulation time necessary to preserve 



adiabaticity as the system grows in size. Recently Raman 
transitions have been driven using a mode-locked high power 
pulsed laser at a wavelength of 355 nm, which is optimum 
for ^^^Yb''" wherein the ratio of differential AC Stark shift to 
Rabi frequency is minimized and spontaneous emission prob- 
abilities per Rabi cycle are < 10^^ per spin [24.1 . 

With this system, it is possible to engineer different Ising 
coupling patterns by controlling the Raman beatnote detun- 
ing /i and observe interesting spin ordering and phase tran- 
sitions some of which can be very sharp, or first order (Esl . 
Long range interactions in this spatially one dimensional sys- 
tem allow for simulating multi-dimensional spin models by 
selectively exciting vibrational modes using multiple Raman 
beatnote detunings. With additional laser beams this scheme 
can potentially simulate more complicated and higher spin- 
dimensional Hamiltonians like XY and XYZ models, which 
map onto nontrivial quantum Hamiltonians such as the Bose- 
Hubbard Hamiltonian. 
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METHODS 
A. Detection of spin states. 

The spin states are detected by spin-dependent fluorescence 
signals collected through //2.1 optics by a photomultiplier 
tube. Spin state | tz) is resonantly excited by the 369.5 
nm detection beam and fluoresces from states, emit- 

ting Poisson distributed photons with mean ~ 12 in 0.8 ms. 
This state appears as 'bright' to PMT. The detection light is far 
off-resonant to spin state | and this state appears 'dark' to 
the PMT. However, due to weak off -resonant excitation bright 
state leaks onto dark state, altering the photon distribution 
ll26ll . Unwanted scattered light from optics and trap electrodes 
also alter the photon distribution. We construct the basis func- 
tion for s bright ions by convolution techniques, and include 
a 5% fluctuation in the intensity of detection beam, which is 
representative of our typical experimental conditions. We then 
fit the experimental data to these basis functions, and obtain 
probabilities P{s). Mean photon counts for dark (m/^) and 
bright states (m^) are used as fitting parameters so as to min- 
imize the error residues. 

The error bars to the order parameters are estimated conser- 
vatively in the following way. Say, best fitting is obtained for 
m£, = ifiD and ms — rhs- The fitting error in those quanti- 
ties are 6m o and <5tob respectively. We generate histograms 
of the order parameters by extracting P{s) ^ 400 times with 
mean dark and bright state counts chosen randomly (Monte 
Carlo method of error analysis) from a Gaussian distribution 
with mean mo, tjib and standard deviations Smn, SniB re- 
spectively. The histogram is fitted to a Gaussian and the width 
is chosen to represent the random error due to fit. The uncer- 
tainty in amount of fluctuation of the detection beam power 
during the experiment is conservatively included in the error 
analysis by repeating the fitting process for a range of fluctu- 
ations. The finite width of the detection beam is taken care of 
by modeling the Gaussian beam having a three step intensity 
profile with appropriate intensity ratios. 



B. Scaled Order Parameters. 

To characterize the spin orders we use different order pa- 
rameters in the experiment, namely the average absolute mag- 
netization per site (rrix) and the Binder cumulant (g). When 
the spins are polarized along the y— direction of the Bloch 
sphere, the distribution of total spin along X— direction is Bi- 
nomial and approaches a Gaussian (with zero mean) in the 
limit of iV ^ cxo. For system size of A^, uix takes on theo- 
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N 

retical value of m° ^ = -^^^ ^ ^C„|Af — 2n| in the per- 

n=0 

feet paramagnetic phase (B/\J\ oo) and unity in the other 
hmit of B/\J\ — 0. In Fig. [Sh and [3}; we rescale mx to 
fhx — {m'^ jY ~ '^k) / w ^ l) which should ideally be 
zero in perfect paramagnetic phase and unity in perfect ferro- 
magnetic phase for any N. This accounts for the 'trivial' finite 
size effect due to the difference between Binomial and Gaus- 
sian distribution. Similarly the Binder Cumulant g is scaled to 
9 = {9% - 9) I {9I - 1) in Fig- where 5^ = 3 - 2/7V is 
the theoretical value of g for J| 00. 



C. Quantum Monte-Carlo Simulations. 

Quantum trajectories are generated by numerically integrat- 
ing the Schrodinger equation, with Hamiltonian (1), while si- 
multaneously executing quantum jumps to account for sponta- 
neous emission and decoherence. Spontaneous emission from 
ion i either localizes the spin of the ion, projecting it into 
2S-yi2\F = 0,mF = 0) (spin state | i,)) or 2Si/2\F = 
IjUiF = 0) (spin state | tz))' or it projects the ion into 
'2Si/2\F — l,TOi? = 1), in which case ion i is factored 
out of the Schrodinger evolution, though it is counted as spin 
up at the time of measurement. Decoherence is modeled 
by the quantum jump operator dx', thus a jump for ion i. 



IV') cr^l'0)> introduces a tt phase shift between the spin 
states I t) and | 4,) (in X— basis) . Jump rates are taken to be 
fixed and equal for all ions. To determine the entangled state 
of the spin ensemble after a spontaneous emission, e.g. from 
ion i, we assume that the ground state configuration prior to 
emission, 

\iz)^'[[h^J\tz)J+S^,J\U)J) 

is mapped, by the far detuned Raman beams, into a very 
small excited-state contribution to the overall system entan- 
gled state, 

A|2Pi/2>. n + 7^j)l tz)j + {^,3 + <5^,j)| U)])], 

with A <C 1 proportional to the amplitude of the Raman beams 
and inversely proportional to their detuning. The (unnormal- 
ized) state after the emission is 

where is | tz)*, | or the factored state 2Sii2\F — 
l,mF = 1). 



